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We have studied experimentally and numerically temperature profiles and the formation of hot 
spots in intrinsic Josephson junction stacks in Bi 2 Sr 2 CaCu2 08 (BSCCO). The superconducting 
stacks are biased in a state where all junctions are resistive. The formation of hot spots in this 
system is shown to arise mainly from the strongly negative temperature coefficient of the c-axis 
resistivity of BSCCO at low temperatures. This leads to situations where the maximum temper- 
ature in the hot spot can be below or above the superconducting transition temperature T c . The 
numerical simulations are in good agreement with the experimental observations. 

PACS numbers: 74.50. +r, 74.72.-h, 85.25.Cp 



I. INTRODUCTION 

Joule heating is an omnipresent issue in current- 
carrying structures and has been studied for a long time. 
General aspects, like the propagation of switching waves 
or the formation of static electrothermal domains in 
bistable conductors are well-known phenomena [HE]- In 
Josephson junctions heating often is small enough to be 
neglected. An exception are stacks of intrinsic Josephson 
junctions (IJJs) in the high temperature superconduc- 
tor Bi 2 Sr 2 CaCu 2 08 (BSCCO). Here, the BSCCO crystal 
structure intrinsically forms stacks of Josephson junc- 
tions, each having a thickness of 1.5 nm. A single IJJ 
may carry a voltage V of some mV and a current / of 
several mA. Although the dissipative power generated by 
a single IJJ is only some /iW, the power inside a stack 
of, say, 1000 IJJs amounts to several mW, with power 




Figure 1: (Color online) Typical design of BSCCO IJJ mesas. 



* Electronic address: hbwangl000@gmail.com 
t Electronic address: kleiner@uni-tuebingen.de 



densities well in excess of 10 4 W/cm 3 . For small sized 
(~ a few jim in diameter, consisting of some 10 IJJs) 
stacks the corresponding overheating has been discussed 
intensively in literature [3HT0], 

Recently, coherent off-chip THz radiation with an ex- 
trapolated output power of some fiW was observed from 
stacks of more than 600 IJJs, with lateral dimensions 
in the 100 /am range [11]. The IJJ stacks have been pat- 
terned in the form of mesa structures, as shown schemat- 
ically in Fig. [T] THz radiation emitted from such IJJ 
stacks became a hot topic in recent years, both in terms 
of experiment [TTH27] and theory [28H56] . 

For these mesas there are two regimes where emis- 
sion occurs p~4j [20] . At moderate input power ( "low-bias 
regime") there is only little heating (< 10 mW), and the 
temperature distribution in the mesa is roughly homoge- 
neous and close to the bath temperature T5. The THz 
emission observed in this regime presumably can be de- 
scribed by more or less standard Josephson physics. At 
high input power ("high-bias regime") a hot spot forms 
inside the mesa p~4j [T9l 120] • The hot spot effectively 
separates the mesa into a "cold" part, which is supercon- 
ducting, and a hot part, which is in the normal state. The 
"cold" part of the mesa is responsible for THz generation 
by the Josephson effect. The hot spot also seems to play 
a role for synchronization [19 , 20, 27 j. It has been found 
that the size and position of the hot spot, and in con- 
sequence also the THz emission, can be manipulated by 
applying proper bias currents across the mesa [19]. Thus, 
in order to understand the mechanism of THz radiation 
in IJJ mesas, it seems crucially important to develop a 
detailed understanding of the hot spot formation. The 
present paper is devoted to this subject. 

In a standard superconducting structure (e.g. a thin 
film) under a strong enough transport current somewhere 
in the sample the resistance rises from zero to a finite 
value, leading to local heating and the formation of a 
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Figure 2: (Color online) (a) Temperature dependence of the 
c-axis resistivity p c , as measured for a 330 x 50 jim 2 wide and 
0.7 /im thick sample for T > T c = 83 K (black circles). For 
lower T, p c has been extrapolated by fitting the IVC, mea- 
sured at Tb — 20 K, using the full 3D heat diffusion equation, 
cf. Sec. IV (b) Temperature dependence of the BSCCO in- 
plane (K a b) and out-of-plane (k c ) thermal conductivity [59] . 



hot spot. To obtain THz emission, IJJ stacks are typ- 
ically biased in a state where all junctions are in their 
resistive state. Here, the out-of-plane resistance R c de- 
creases continuously when heating the sample through 
T c [57l|58], cf. Fig. [5] (a). In-plane currents still flow 
with zero resistance below T c and with finite resistance 
above T c . However, even in the normal state these layers 
add only a minor contribution to the total voltage across 
the IJJ stack and thus to the overall power dissipation 
due to the huge ratio p c /p a b > 10 5 of the out-of-plane to 
the in-plane resistivity. It is unlikely that this contribu- 
tion gives rise to hot spot formation. Also the BSCCO 
thermal conductance varies relatively weakly with tem- 
perature [59], cf. Fig. [2] (b) . Thus, the above mechanism 
of hot spot generation does not work and T c is no longer 
a peculiar temperature for the thermal balance of the 
sample. 

There are other ways to create hot domains in sys- 
tems which may or may not be related to superconduc- 
tivity [TJ[2]. In particular, current- voltage characteristics 
(IVCs) and the thermal breakdown were studied in sys- 
tems having a resistivity decreasing with increasing tem- 
perature (negative-temperature-coemcient resistor) [60] - 
[62] . The IVCs of these resistors strongly resemble the 
IVCs measured for IJJ mesas. Especially, the appearance 
of a hot domain leads to an abrupt change in differential 
resistance. The quantity in common, a strongly negative 
dR/dT, is the key to understand hot spot formation in 
BSCCO mesas. 

Recently, Yurgens et al. have simulated the thermal 
heating and the temperature distribution in BSCCO IJJ 
mesas [47, 49 , using a 3D finite-element software [63]. In 
this pioneering work the electrical and thermal proper- 
ties of the various current carrying and insulating layers 
were taken into account. The formation of hot spots ob- 
served in [14] was reproduced qualitatively. However, the 
occurring phenomena need further study. For example, 
the IVCs in [47] [49] have been calculated using a self- 
consistent procedure based on Newton's law of cooling 
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Figure 3: (Color online) Discrete approximation for a mesa, 
(a) Dimensions of mesa and base crystal, (b) The mesa is re- 
placed by two vertically cooled resistors Ra and Rb producing 
Joule heat Qa and Qb, which is vertically transported to a 
thermal bath via heat-transfer powers Wa and Wb> 



and Ohm's law and do not exhibit the experimentally 
observed abrupt changes in differential resistance when 
the hot spot appears. They resemble much less the ex- 
perimental curves than the ones calculated in [60H62] . 

A complete study of the Josephson effect in BSCCO 
mesas in the presence of hot spots is a formidable and un- 
solved issue. In this paper we are treating experimentally 
and theoretically hot spot formation in BSCCO mesas. 
In the theoretical part of our study the presence of the 
Josephson effect, i.e. THz radiation, the formation of 
electromagnetic standing waves, interactions between hot 
spots and waves etc. is not considered. This approach 
to hot spot formation seems justified, since the emitted 
radiation power is 3-4 orders of magnitude lower than 
the dc input power. It may, however, serve as a zero or- 
der approximation towards solving the full problem. In 
the simulations we derive the electrical current density in 
the mesa under investigation and thus also the potential 
difference between top and bottom electrodes, directly 
generating the IVC for a sample, following [60-62 rather 
than gang. 

The paper is organized as follows. In Sec. [H]we con- 
sider a simple discrete resistor model to get a basic under- 
standing of the heating phenomena involved. In Sec. Ill 
a ID model is discussed which is extended to 3D and real- 



istic sample geometries in Section IV The discussions in 
these sections are based on the thermal and electrical pa- 
rameters of the BSCCO crystals, as used in experiment. 
In Sec. |lV|we also address experimental observations, as 
made in~Pl [19] [20] [27] . Sec. |V| concludes our work. 



II. DISCRETE RESISTORS 

The electrothermal behavior of conducting materials 
can be investigated by considering the heat balance equa- 
tion between Joule self- heating Q(T, A) and the heat 
transfer power W(T) to the coolant Q(T, A) = W(T) pQ. 
Here A is some control parameter (in our case the volt- 
age V across the sample). To approach the experimental 
situation of IJJ mesas we first briefly study the model 
of two current-biased resistors Ra(Ta) and Rb(Ta) con- 
nected in parallel, each representing one half of a mesa 
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of length Z, width w and height h, cf. Fig. [3j Ta and 
Tb are the temperatures of these resistors. Ra and Rb 
shall be equal for T4 = Tb- Joule heating is produced 
via Qi — IiVi, where i = (A,B). The total current 
is / = I a + Ib and further Va = Vb, i.e. we neglect 
the voltage drop due to in-plane currents. The resis- 
tors are thermally connected to a bath (temperature X5), 
which, at a distance L (the thickness of the base crys- 
tal) removes heats Wa and Wb "vertically" , through the 
BSCCO out-of-plane thermal conductivity k c . We first 
assume Ta — Tb — T. Then, the IVC of the mesa can 
be parameterized by T, using Q = W [64] : 



V = y/R(T)W(T-T h ); I 



> W{T-T b ) 
R(T) 



(1) 



with W(T - T b ) = (lwK c /L) • (T - T b ) and R(T) = 
(h/lw) • p c (T). For further calculations we use a con- 
stant k c = 0.6 W/mK. As we want to study the question 
whether or not the particular p c (T) of our mesas can lead 
to hot spots we use a temperature dependence which is 
as close as possible to the experimental situation. Above 
T c we obtained p c (T) from the out-of-plane resistance of 
one of our mesas, cf. solid circles in Fig. [2] (a). Be- 
low T c , p c (T) is extrapolated by fitting the measured 
IVC of the mesa at a bath temperature of 20 K (see be- 
low), using the full 3D heat diffusion equation (dashed 
line in Fig. [2] (a)). L = 17 pm is chosen, which is a 
typical value for the thickness of the BSCCO base crys- 
tal of the samples we want to discuss p~4j [19j [20j [27] . 
Length, width and height of the mesa are, respectively, 
taken to be 330, 50 and 1 /im, representing sample 1 from 
[20] . With these dependencies, the calculated IVC of the 
mesa is S-shaped and shows a region of negative differ- 
ential resistance, cf. solid line in Fig. [4] (c). In this volt- 
age region thermal bistability can occur, since W = Q 
holds for more than one value of T [1 . In fact, writ- 
ing dV/dl = (dV/dT)/(dI/dT) < 0, using Eq. and 
W oc (T - T h ) one obtains — (T - T h ) (dR/dT)/R > 1 as 
a condition for obtaining negative differential resistance 
in the IVC and thus the possibility to have thermal in- 
stability [2 . 

We assume for the following, that / and thus Q is 
increased from zero step-by-step. Fig. [4] (a) and (b) 
show the individual IVCs of resistor A and B respectively, 
while (c) shows the IVC of the whole mesa. For small Q 
the temperature is the same in both resistors and they 
carry the same current. In principle, further increase of 
I would make the whole mesa pass the point S of local 
maximal voltage Vb, cf. Fig. [4] (c) and enter the unstable 
PQ area of negative differential resistance. This is exem- 
plarily indicated for point a in Fig. [4](c). Here, the two 
resistors with equal temperature T a , would be in states 
a a and cf. Fig. [4] (a) and (b) respectively. The insta- 
bility and the constraints of equal voltage and fixed cur- 
rent force the mesa into the state /?, which is composed of 
state Pa with T$ A > T a and /3b with T$ B < T a , cf. Fig. 
[4] (a) and (b) respectively. The combination of /3a and 
/3b is the only stable solution. The resulting total IVC of 
Fig. [4] (c) follows the path indicated by the dashed (red) 
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Figure 4: (Color online) Hot spot formation in a two-resistor 
model, cf. Fig. |3jb). (a) and (b) display the IVCs of the 
two individual parts A and B, respectively, (c) shows the 
IVC of the combined system. The axes are normalized to the 
current (voltage) of the point showing local maximal voltage 
Vb- The total current through the mesa at Vb is Jo- The 
bias points indicated by Greek characters are discussed in the 
text. In (c) for the solid (blue) curve resistors A and B are at 
the same temperature, while for the dashed (red) curve their 
temperature differs, corresponding to hot spot formation in 
the continuous case. 



line, differing in voltage from the isothermal case (solid 
blue line). With increasing /, starting from point J, the 
points Pa and /3b "move" towards lower voltage. Note 
that this implies, that the cold resistor becomes colder 
while the hot resistor keeps increasing its temperature. 
When /3a has reached the minimal voltage, both /3a and 
(3b start to move towards larger voltage, i.e. also the 
temperature of the cold part starts to increase. Finally, 
when /3b reaches the voltage Vq, Ta ^ Tb becomes im- 
possible and the mesa switches back to the homogeneous 
solution. 

The model of two parallel resistors can be extended by 
several ways: First, an in-plane thermal coupling Wab 
between resistors A and B may be included. Then, the 
cold part will cool the hot part and (thermal) differences 
between A and B will be less severe. This will shift the 
point, where the homogeneous solution and the solution 
Ta 7^ Tb fork, to higher input power [60, 62 . Also, the 
difference in voltage between the homogeneous and the 
inhomogeneous solution will be diminished [60j [62]. A 
detailed discussion, however, is out of the scope of this 
section. In-plane cooling will be taken into account in 
the subsequent sections. Second, one may allow the two 
resistors (the area of the "hot" and "cold" parts) to be 
unequal and variable in size. Then, one faces a continu- 
ous set of solutions. Third, one may consider more than 
two resistors in parallel. This would be also applicable to 
the description of arrays of IJJ stacks, which are inter- 
esting for obtaining a large THz emission output power. 
In this scenario, the whole system will tend to a state, 
where only one of the stacks is hot, while all the others 
are cold 62 . 
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III. ID MODEL 

In this section we consider a ID continuous model to 
find the temperature distribution in the mesa for the 
simplest continuous case, still treating hot spot forma- 
tion from a generic point of view. That is, we assume 
a thin (along z) and narrow (along y) mesa, neglecting 
T-variations along z- and ^-direction in the mesa (see 
[TJ [62] for details). Then, T = T(x) is defined by the 
heat diffusion equation: 



-h- 



d 

dx 



«c (T) 



(T - T b ) 



V 



Pc (T) h 



(2) 



The first term describes the thermal diffusion in x- 
direction and the second one the cooling due to the base 
crystal with the coefficient K c jL regulating its strength. 
The third term represents Joule heating. The sample di- 
mensions L, /i, I and w are defined in Fig. [3] (a). We use 
L = 19 //m, h = 1 /mi, / = 330 //m, w = 50 //m and ^,5, n c 
and p c as in Figs. [2] (a), (b). The boundary conditions are 
chosen to be dT/dx (x = 0) = dT /dx (x = I) = 0. These 
boundary conditions neglect edge cooling. To solve Eq. 
p| for a given current / we use V = Ih/ J p~ x (T) dxdy. 
We numerically solve Eq. ([2| using finite element analysis 
[63] . Note that there is always a homogeneous solution. 
To find a nontrivial T (x), a proper initial function Ti (x) 
has to be used. A calculated IVC for T5 = 20 K is shown 
in Fig. [5] (a). It resembles the shape of the IVC of the 
two-resistor model, cf. Fig. [4] Figure [5] (b) shows the 
temperature in the mesa for the homogeneous solution 
at the bias points indicated in (a). One notes that the 
mesa temperature is below T c up to quite high currents 
~ 40 mA. Figure [5] (c) shows solutions for the bias points 
indicated in the IVC, when a hot spot has formed. Here, 
the temperature in the hot part rises rapidly to temper- 
atures well above T c , while the temperature of the cold 
part is near T5 = 20 K. Also, one observes, that the hot 
part grows in size when / is increased. Further note, that 
in the presence of a hot spot the temperature T co u of the 
cold part is below the temperature of the homogeneous 
solution for the same value of Q; T co id decreases with in- 
creasing Q, and finally converges against a limiting value. 
The strength of the deviation of the temperature profile 
from the homogeneous solution directly correlates with 
the strength of branching in the IVC. In the depicted 
case the branching is very strong, which is due to a small 
ratio of the in-plane to the out-of-plane thermal coupling, 
cf. first and second term in Eq. 

For a given / the hot spots presented in Fig. [5] (c) are 
not the only possible solutions to Eq. Q [65 . For sym- 
metry reasons also the mirrored solution exists, as well as 
solutions with the hot spot near and in the center of the 
mesa, cf. Fig. [5](d). In the IVC the different solutions 
slightly differ in V and can be traced over some range in /. 
Thus, the IVC consists of several branches distinguishing 
specific kinds of hot domains. Experimentally, in some 
cases, hot spot formation in different places of the mesa 
has been detected by low-temperature scanning-laser mi- 
croscopy (LTSLM). However, usually a specific configu- 
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Figure 5: (Color online) Simulation results of Eq. |2| for 
Tb — 20 K and L — 19 /mi. Other mesa dimensions are listed 
in Sec. [HI Graph (a) shows the IVC for the homogeneous 
solution (black curve) and a solution showing a hot spot on 
the right mesa end (gray curve). T(x)-profiles are displayed 
in (b) for the homogeneous case and in (c) for the hot spot 
case. The numbers indicate the bias points on the IVC. In 
(d) T(x)-profiles, obtained from different Ti(x), are shown 
for solutions with / = 9.5 mA, exhibiting various shape and 
positioning of the hot spot. 



ration is much more stable than the others, presumably 
due to inhomogeneities like attached wires. In the cal- 
culations, also solutions with more than one hot domain 
[66] can be found, cf. Fig. |5](d). However, this has not 
been observed in any of our LTSLM measurements. It is 
argued in [62 , that such a state is very unstable and will 
not occur, since the sample can be seen as a parallel cir- 
cuit of several discrete parts with small thermal coupling 
between them, cf. section [TTl 



IV. 3D MODEL 

In this section we address hot spot formation in 3D. 
The goal is to quantitatively compare our experimental 
observations with the numerical simulations, using the 
same code [63] as in [47| |49] . Similarly, we also include 
various electrically and thermally conducting and insu- 
lating layers that are in contact with our BSCCO mesas. 
The electrical, thermal and geometrical parameters used 
for the calculations are as close as possible to the exper- 
imental situation. The geometry used is still somewhat 
simplified compared to the real samples but should al- 
low to capture the relevant physics. Figure [6] depicts the 
model. The substrate is omitted and a boundary condi- 
tion T5 = const, is applied to the bottom surface of the 
glue layer, representing the thermal bath. This simplifi- 
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Figure 6: (Color online) Model geometry of the mesa. 



cation can be done with very little impact on the results 
for the mesa, since the thermal conductivity of the sub- 
strate (e.g sapphire) is by far better than that of the glue 
layer. The geometric dimensions of the mesa, the thick- 
nesses of the glue layers (10 - 30 /im) and of the gold 
coatings {Hau ~ 30 nm) were roughly chosen as in the 
real samples. The base crystal's lateral size is typically of 
the order of 1 mm, while its thickness hb c may vary from 
about ten to several hundred /im, strongly depending on 
the fabrication process. The current leads are simply rep- 
resented by boundary conditions on the surfaces of either 
the gold layer on the mesa or on the base crystal. The 
current / is injected through a 20 x 10 /im 2 rectangle 
and the current sink is defined as a ground of large area 
(roughly 0.3 mm ), cf. Fig. [§ The voltage across the 
mesa is obtained as the potential difference between the 
two electrodes. 

The equation to be solved is [49] : 

-V[«(T(r)) VT(r)]=p(r(r)) j 2 (r) , (3) 

where p and n are the resistivity and thermal conductiv- 
ity tensor, respectively, and r is the spatial coordinate. 
Unlike the mesa, the base crystal is not always in the 
resistive state. We model its resistance by using the p c 
vs. T data indicated by solid circles in Fig. [5] (b). The 
in-plane resistivity p ab is the same for both mesa and 
base crystal; we use the same T dependence as in [49] . 
The thermal conductivity for BSCCO is used from [59] , 
cf. Fig. [5] (a). Thermal and electrical conductivity for a 
30 nm thick Au film are adopted from [67] . For the ther- 
mal conductivity n g \ ue of the glue between the BSCCO 
base crystal and the substrate to first order we use the 
polyimide data of [68 . Since our glue might have slightly 
different properties we in addition multiply K g \ ue with a 
factor n g i ue , which we fit by adjusting the calculated IVC 
to the measured one. 

The base crystal introduces an effective side-cooling 
of the mesa, which in general makes a solution showing 
variation in x and y direction (with or without hot spots) 
favorable. Indeed, in contrast to the one-dimensional cal- 
culations, hot spot solutions appeared basically by them- 
selves, i.e. it was not necessary to find them by choos- 
ing a proper initial condition. The side-cooling leads to 
an elliptic shape of the hot spot (for rectangular shaped 
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Figure 7: (Color online) Comparison of 3D simulation and ex- 
perimental data for sample 1 from [20] at T& — 20 K. (a) shows 
the measured (black, solid circles) and simulated (red solid 
line) IVCs. In (b) simulated T(x)-profiles along the dashed 
line indicated in Fig. [6]at z — 0.5 h are shown. The diamonds 
in (a) indicate the corresponding bias points. The calculated 
and measured AV(x) are shown in (c) and (d), respectively. 
Diamonds in (c) indicate the x-position where T — T c . 



mesas). Also, the hot spot is not limited to the mesa 
itself anymore, but may extend significantly in lateral di- 
rection into the base crystal (see below). This is exactly 
what has been found experimentally [19]. The same oc- 
curs in z-direction, as has been discussed in [49] . 

We investigate sample 1 from [20]. The electrical and 
thermal parameters of this sample have already been used 
in the previous sections. We have further used the pa- 
rameters: h bc =40 /im, h g i ue — 25 /im, lgi ue — 1 mm and 
rigiue = 1-95. Figure [7] (a) compares the measured IVC 
with the calculated one for T5 = 20 K. The good agree- 
ment stems from the fact that we have adjusted <j c (T) 
below T c and n g i ue to match this curve. The simulated 
T(x)-profiles, calculated along the dashed line in Fig. [6] 
at z = 0.5 /i, are shown in Fig. [7] (b) for the bias points 
indicated in Fig. [7] (a). They show an almost constant 
temperature in the low bias regime, whereas for increas- 
ing current the hot spot forms by the growth of a buckling 
in the T(x)-profile (compare curves 2 and 3 in Fig. 0(b)). 
Further increasing / and thus Q leads to a growth in di- 
ameter and maximal temperature of the hot spot (curves 
4 to 8). Note that T c , indicated by the horizontal dashed 
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line, can be significantly exceeded in the center of the hot 
spots, confirming the results in [49] . 

We next want to provide a quantitative comparison be- 
tween the hot spot signals observed in LTSLM [T4l[T9ll2Q] 
and the calculated temperature distributions for this 
sample. In LTSLM the laser spot at position (xo,yo) 
causes a maximum temperature rise AT ~ 1-3 K, de- 
pending on the laser power. In turn there is a change 
AV(xo,yo) in the voltage V across a sample. One often 
has a response which partially arises from the reduction 
of the Josephson critical current density and partially 
from the change in resistance, see e.g. [69 j. However, if 
dR c (T)/dT dominates the thermal physics, AV(xo,yo) 
can be treated as in [70] yielding 



AV(x 0l y ) 



-IR 2 eff ATA L da c 



h 



dT 



(T(x ,y ))- (4) 



R e ff = V/I is the (ohmic) sample resistance at a 
given I and Al is the effective area warmed up by the 
laser (some /im 2 ). dcr c /dT (T (xq,2/o)) denotes the tem- 
perature derivative of the c-axis electrical conductivity. 
The calculated and measured AV taken at various bias 
points indicated in Fig. [7] (a), are shown in Fig. [7] (c) 
and (d), respectively. For the simulations we have used 
AT • Al = 56K/im 2 . The value makes sense, since we 
expect a temperature rise AT ~ 2K and Al ~ 25 /im 2 
for the samples we discuss here. The calculated curves 
agree reasonably well with the measurements, although 
differences occur at low bias and near the hot spot nu- 
cleation point. Particularly, for the bias points 1 and 2, 
the simulation yields a parabolic shape of AV, while the 
experimental data are shaped less regular. Note, how- 
ever, that in these regions the Josephson currents, which 
are neglected in our analysis, may play a major role. For 
curve 3 in the simulation hot spot formation has already 
occurred, while in experiment the mesa is close to the 
nucleation point but still under critical. For a bias well 
above the hot spot nucleation point theoretical curves 
and experimental data agree well. Specifically the dou- 
ble hump feature in A V (x) is reproduced correctly in 
the simulations. The local temperature at the maxima 
in AV corresponds to the temperature T* « 80 K, for 
which d(j c /dT is maximum, cf. Eq. Q. Between the 
two AV maxima, T > T*. By coincidence, T* w T c ; the 
diamonds in Fig. [7] (c) indicate the locations for which 
T = T c . Thus, the border between superconducting and 
non-superconducting parts, which is important for THz 
emission, can be approximately identified by the position 
of the humps. 

We next investigate the dependence of hot spot for- 
mation on T b . Figure [8] shows a similar set of data as 
Fig. [7| but for T b = 42 K. Here, n g i ue = 3.5 has been 
chosen. The transition region between the hot and cold 
domain is less steep than for T b = 20 K. Also, the nucle- 
ation point of the hot spot has moved to higher currents 
(10 mA for T b = 20 K and 14 mA for T b = 42 K) and the 
back-bending of the IVC has decreased. These effects 
arise from the fact that the xy-pl&rie thermal coupling 
has increased relative to the out-of-plane thermal cou- 
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Figure 8: (Color online) Comparison of 3D simulation and ex- 
perimental data for sample 1 from [20] at Tb — 42 K. (a) shows 
the measured (black, solid circles) and simulated (red solid 
curve) IVC. In (b) simulated T(x)-profiles along the dashed 
line indicated in Fig. [6]at z — 0.5 h are shown. The diamonds 
in (a) indicate the corresponding bias points. The calculated 
and measured AV(x) are shown in (c) and (d), respectively. 
Diamonds in (c) indicate the x-position where T — T c . 



pling 



cf. Sec. Ill 



Note, that this also means for 
equal input power, that the hot domain reaches higher 
and the cold domain reaches lower temperatures for T b = 
20 K as for T b = 42 K. Figures [8] (c) and (d) respectively 
show the calculated and measured LTSLM profiles. As 
for [7] (c) and (d) the agreement between experimental 
data and simulations is reasonable, except for the bias 
point where hot spot formation sets in (curve 4) . For the 
calculations we have used AT • Al = 16K/im 2 , which is 
by a factor of 3 lower than for the case of T b = 20 K. This 
is attributed to a reduced incident laser power, which had 
been readjusted for every measurement. 

In Fig. [9] (a), we show the T (x) -profile, calculated 
along the dashed line in Fig. [6] at z = 0.5 h, for 3 values 
of T b . For all curves, V = 0.8 V. This condition has been 
motivated by measurements of the linewidth Af of THz 
radiation [27]. Here, for Af vs. T5, taken at a fixed emis- 
sion frequency (corresponding to V = const, for a fixed 
number of oscillating IJJs) the dependence Af oc T 6 -4 
has been found unexpectedly. We are interested in the 
question whether or not corresponding changes with T b 
can be seen in the T distribution in the mesa. Figure[9](a) 
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x (jiim) / (mA) 

Figure 9: (Color online) (a) Simulated T(x)-profiles along the 
dashed line indicated in Fig. |6]at z = 0.5 h for three different 
values of Tb at constant V — 0.8 V. The dc power P — IV is 
26.4, 18.2 and 14.4 mW for T b = 20, 30 and 42 K, respectively, 
(b) P vs. / at 20 K for measurement (green squares), sim- 
ulation with homogeneous T (red dashed line) and hot spot 
formation (blue, solid line). The pink circles depict experi- 
mental data and the black, short-dashed line simulated values 
with hot spot formation at 42 K. 




Figure 10: (Color online) (a)-(d) Surface plot of hot spot 
solutions, obtained by Eq. ^ for a mesa with two current 
injection points, indicated by black rectangles. The sum Io 
of the currents through the left (Ii) and right (I r ) injection 
points has been kept constant, and for the ratio Ii/Iq values 
of (a) 1, (b) 0.7, (c) 0.5 and (d) 0.425 have been used. Graph 
(e) shows the center position Xh of the hot spot vs. current 
injection ratio for simulated and measured data. 



shows, that the peak temperature in the mesa is higher 
at low T5 than at high while the coldest temperatures, 
reached at the right edge of the mesa, behave oppositely. 
Thus, thermal gradients at low X5 are stronger than at 
high T5. However, this effect roughly changes linearly 
with T5 and presumably cannot explain the Af ex T 5 -4 
dependence. 

Figure [9] (b) compares for two values of T5 the mea- 
sured and simulated DC power P = IV as a function of /. 
One observes two regimes, each with a roughly constant 
slope. The first - low-bias - regime has no hot spot and, 
for T b = 20 K, spans from to 10 mA (14 mA for 42 K), 
whereas the second - high-bias - regime has a hot spot 
and begins at 10 mA (14 mA for 42 K). Interestingly, 
at the intersection of these two regimes the maximum 
temperature in the mesa has reached the temperature 
fulfilling da c /dT = 0. This point also corresponds to the 
kink in the IVC, observed for several mesas. Note that 
the calculation for homogeneous T (red dashed curve) 
shows no such kink. A plot like this may thus be helpful 
to distinguish in an experiment, whether or not one has 
reached the regime with hot spots. 

The last issue we want to address is the correlation 
between the point of current injection and the location, 
where the hot spot is established. Typically, in exper- 
iment the appearance of the hot spot was close to, but 
not exactly at the bond wire to the mesa surface [19 j. We 
see the same effect in our simulations, Fig. 10 (a) illus- 
trates this for a situation, where the current is injected 
from the left. Here, the side-cooling prevents the hot spot 
from nucleating at the very left end of the mesa, result- 
ing in a positioning of the hot spot at several fim right of 
the current injection point. Further, it has been shown 
that by using two injectors located on opposite sides of 
the mesas the hot spot can be moved by changing the 



ratio of currents through these injectors [19]. Figure 10 



(a)-(d) show a sequence of calculations where the ratio 



between injection currents through the left (current I\) 
and right (current I r ) was varied, keeping the sum of 
the currents Iq constant. We used ratios Il/Io of, re- 
spectively, 1, 0.7, 0.5 and 0.425. As one can see, the 
hot spot indeed can be moved continuously, as in experi- 
ment. In Fig. 10 (e) we have plotted the center position 
of the hot spot as a function of Ii/Iq. Experimental data 
are shown by (red) squares and theoretical data by the 
(black) solid line. The agreement is reasonable, showing 
that this effect can be essentially understood from the 
thermal calculations presented in this paper. 

Finally, we briefly mention that also two other geome- 
tries discussed in [19 can be reproduced very well in the 
3D simulation - a disk shaped mesa and a mesa of Y 
shape, where the hot spot forms at the intersection of 
the three lines, although the bias current injection point 
was at the foot of the Y. 



V. CONCLUSION 

In conclusion, we have investigated experimentally and 
numerically the temperature profiles and hot spot forma- 
tion in IJJ mesas. We have shown, that the hot spots 
dominantly arise because of the strongly negative tem- 
perature coefficient of the out-of-plane resistance of the 
mesas. This mechanism is different from the more con- 
ventional hot spot formation in superconductors and, in 
particular, allows for hot spots with a maximum temper- 
ature below as well as above the transition temperature 
T c . We have given - in the frame of what available data 
allow - a quantitative comparison between simulation 
and experiment, showing reasonable agreement. Numer- 
ous effects observed in previous papers on hot spot for- 
mation in intrinsic Josephson junction stacks [14} [19} [20] 
are reproduced by the simulations, making us confident 
that the description given in this paper captures the es- 
sential physics, except for the interplay of hot spots and 
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THz waves. Resolving this issue is a task for the future. 



Acknowledgments 



operative Program and the German Israeli Foundation 
(Grant No. G-967-126.14/2007). S.G. acknowledges 
support by AFOSR. 



We gratefully acknowledge financial support by the 
JST/DFG strategic Japanese-German International Co- 



[1] A. V. Gurevich and R. G. Mints, Rev. Mod. Phys. 59, 
941 (1987). 

[2] A. F. Volkov and S. M. Kogan, Sov. Phys. Uspekhi 11, 
881 (1969). 

[3] P. Thomas, J. Fenton, G. Yang, and C. Gough, Physica 

C 341, 1547 (2000). 
[4] J. C. Fenton and C. E. Gough, J. Appl. Phys. 94, 4665 

(2003). 

[5] K. Anagawa, Y. Yamada, T. Shibauchi, M. Suzuki, and 

T. Watanabe, Appl. Phys. Lett 83, 2381 (2003). 
[6] V. N. Zavaritsky, Phys. Rev. Lett. 92, 259701 (2004). 
[7] A. Yurgens, D. Winkler, T. Claeson, S. Ono, and 

Y. Ando, Phys. Rev. Lett. 92, 259702 (2004). 
[8] V. M. Krasnov, M. Sandberg, and I. Zogaj, Phys. Rev. 

Lett. 94, 077003 (2005). 
[9] H. B. Wang, T. Hatano, T. Yamashita, P. H. Wu, and 

P. Miiller, Appl. Phys. Lett. 86, 023504 (2005). 
[10] B. Verreet, N. Sergeant, D. M. Negrete, M. Torstensson, 

D. Winkler, and A. Yurgens, Supercond. Sci. Technol 20, 

S48 (2007). 

[11] L. Ozyuzer, A. E. Koshelev, C. Kurter, N. Gopalsami, 
Q. Li, M. Tachiki, K. Kadowaki, T. Yamamoto, H. Mi- 
nami, H. Yamaguchi, et al., Science 318, 1291 (2007). 

[12] K. Kadowaki, H. Yamaguchi, K. Kawamata, T. Ya- 
mamoto, H. Minami, I. Kakeya, U. Welp, L. Ozyuzer, 
A. Koshelev, C. Kurter, et al., Physica C 468, 634 (2008). 

[13] L. Ozyuzer, Y. Simsek, H. Koseoglu, F. Turkoglu, 
C. Kurter, U. Welp, A. E. Koshelev, K. E. Gray, W. K. 
Kwok, T. Yamamoto, et al., Supercond. Sci. Technol. 22, 
114009 (2009). 

[14] H. B. Wang, S. Guenon, J. Yuan, A. Iishi, S. Arisawa, 
T. Hatano, T. Yamashita, D. Koelle, and R. Kleiner, 
Phys. Rev. Lett. 102, 017006 (2009). 

[15] H. Minami, I. Kakeya, H. Yamaguchi, T. Yamamoto, 
and K. Kadowaki, Applied Physics Letters 95, 232511 
(pages 3) (2009). 

[16] C. Kurter, K. E. Gray, J. F. Zasadzinski, L. Ozyuzer, 
A. E. Koshelev, Q. Li, T. Yamamoto, K. Kadowaki, W.- 
K. Kwok, M. Tachiki, et al., IEEE Trans. Appl. Super- 
cond. 19, 428 (2009). 

[17] K. E. Gray, L. Ozyuzer, A. K. C. K. Kurter, K. Kad- 
owaki, T. Yamamoto, H. Minami, H. Yamaguchi, M. T. 
M, W. Kwok, and U. Welp, IEEE Trans Appl. Super- 
cond. 19, 3755 (2009). 

[18] K. Kadowaki, M. Tsujimoto, K. Yamaki, T. Yamamoto, 
T. Kashiwagi, H. Minami, M. Tachiki, and R. A. Klemm, 
J. Phys. Soc. Jpn 79, 023703 (2010). 

[19] S. Guenon, M. Griinzweig, B. Gross, J. Yuan, Z. Jiang, 
Y. Zhong, A. Iishi, P. Wu, T. Hatano, D. Koelle, et al., 
Phys. Rev. B 82, 214506 (2010). 

[20] H. B. Wang, S. Guenon, B. Gross, J. Yuan, Z. G. 
Jiang, Y. Y. Zhong, M. Gruenzweig, A. Iishi, P. H. Wu, 
T. Hatano, et al., Phys. Rev. Lett. 105, 057002 (2010). 



[21] M. Tsujimoto, K. Yamaki, K. Deguchi, T. Yamamoto, 
T. Kashiwagi, H. Minami, M. Tachiki, K. Kadowaki, and 
R. A. Klemm, Phys. Rev. Lett. 105, 037005 (2010). 

[22] H. Koseoglu, F. Turkoglu, Y. Simsek, and L. Ozyuzer, J. 
Supercond. Nov. Magn. 24, 1083 (2011). 

[23] T. M. Benseman, A. E. Koshelev, K. E. Gray, W.-K. 
Kwok, U. Welp, K. Kadowaki, M. Tachiki, and T. Ya- 
mamoto, Phys. Rev. B 84, 064523 (2011). 

[24] K. Yamaki, M. Tsujimoto, T. Yamamoto, A. Furukawa, 
T. Kashiwagi, H. Minami, and K. Kadowaki, Opt. Ex- 
press 19, 3193 (2011). 

[25] T. Kashiwagi, M. Tsujimoto, T. Yamamoto, H. Mi- 
nami, K. Yamaki, K. Delfanzari, K. Deguchi, N. Orita, 
T. Koike, R. Nakayama, et al., J. J. Appl. Phys. 51, 
010113 (2012). 

[26] M. Tsujimoto, T. Yamamoto, K. Delfanazari, 
R. Nakayama, T. Kitamura, M. Sawamura, T. Kashi- 
wagi, H. Minami, M. Tachiki, K. Kadowaki, et al., Phys. 
Rev. Lett. 108, 107006 (2012). 

[27] M. Li, J. Yuan, N. Kinev, J. Li, B. Gross, S. Guenon, 
A. Ishii, K. Hirata, T. Hatano, D. Koelle , et al., Phys. 
Rev. B 8 6, 06050 5 (2012), URL |http: //link. aps.org/ 
doi/10 . 1 103/PhysRevB . 86 . 060505| 

[28] L. N. Bulaevskii and A. E. Koshelev, Phys. Rev. Lett. 
99, 057002 (2007). 

[29] A. E. Koshelev and L. N. Bulaevskii, Phys. Rev. B 77, 
014530 (2008). 

[30] A. E. Koshelev, Phys. Rev. B 78, 174509 (2008). 

[31] S. Lin and X. Hu, Phys. Rev. Lett. 100, 247006 (2008). 

[32] V. M. Krasnov, Phys. Rev. Lett. 103, 227002 (2009). 

[33] X. Hu and S. Z. Lin, Phys. Rev. B 80, 064516 (2009). 

[34] Y. Nonomura, Phys. Rev. B 80, 140506 (2009). 

[35] M. Tachiki, S. Fukuya, and T. Koyama, Phys. Rev. Lett. 
102, 127002 (2009). 

[36] T. Koyama, H. Matsumoto, M. Machida, and K. Kad- 
owaki, Phys. Rev. B 79, 104522 (2009). 

[37] N. Pedersen and S. Madsen, IEEE Trans Appl. Super- 
cond. 19, 726 (2009). 

[38] R. A. Klemm and K. Kadowaki, Journal of Physics: Con- 
densed Matter 22, 375701 (2010). 

[39] V. M. Krasnov, Phys. Rev. B 82, 134524 (2010). 

[40] S. O. Katterwe, A. Rydh, H. Motzkau, A. B. Kulakov, 
and V. M. Krasnov, Phys. Rev. B 82, 024517 (2010). 

[41] S. Savel'ev, V. A. Yampol'skii, A. L. Rakhmanov, and 
F. Nori, Rep. Prog. Phys. 73, 026501 (2010). 

[42] X. Hu and S. Z. Lin, Supercond. Sci. Technol. 23, 053001 
(2010). 

[43] T. Tachiki and T. Uchida, J. Appl. Phys. 107, 103920 
(2010). 

[44] S. Z. Lin and X. A. Hu, Phys. Rev. B 82, 020504 (2010). 
[45] W. Zhou, C. Wang, and Q.-H. Chen, Phys. Rev. B 82, 

184514 (2010). 
[46] A. E. Koshelev, Phys. Rev. B 82, 174512 (2010). 



9 



[47] A. Yurgens, M. Torstensson, and D. Winkler, Physica C: 

Superconductivity 470, 818 (2010). 
[48] M. Tachiki, K. Ivanovic, K. Kadowaki, and T. Koyama, 

Phys. Rev. B 83, 014508 (2011). 
[49] A. A. Yurgens, Phys. Rev. B 83, 184501 (2011). 
[50] A. A. Yurgens and L. N. Bulaevskii, Supercond. Sci. 

Technol 24, 015003 (2011). 
[51] V. M. Krasnov, Phys. Rev. B 83, 174517 (2011). 
[52] S. O. Katterwe and V. M. Krasnov, Phys. Rev. B 84, 

214519 (2011). 

[53] T. Koyama, H. Matsumoto, M. Machida, and Y. Ota, 
Supercond. Sci. Technol. 24, 085007 (2011). 

[54] S. Z. Lin and X. A. Hu, J. Nanoscience and Nanotechnol. 
11, 2916 (2011). 

[55] S.-Z. Lin, X. Hu, and L. Bulaevskii, Phys. Rev. B 84, 
104501 (2011). 

[56] T. M. Slipchenko, D. V. Kadygrob, D. Bogdanis, V. A. 
Yampol'skii, and A. A. Krokhin, Phys. Rev. B 84, 224512 
(2011). 

[57] A. Yurgens, D. Winkler, N. V. Zavaritsky, and T. Clae- 

son, Phys. Rev. Lett. 79, 5122 (1997). 
[58] Y. I. Latyshev, T. Yamashita, L. N. Bulaevskii, M. J. 

Graf, A. V. Balatsky, and M. P. Maley, Phys. Rev. Lett. 

82, 5345 (1999). 
[59] M. F. Crommie and A. Zettl, Phys. Rev. B 43, 408 

(1991). 

[60] E. Spenke, Electrical Engineering (Archiv fuer Elek- 

trotechnik) 30, 728 (1935). 
[61] H. Lueder and E. Spenke, Physikalische Zeitschrift 36, 

767 (1936). 

[62] E. Spenke, Wissenschaftliche VeroefTentlichungen aus 



den Siemens- Werken 15, 92 (1936). 

[63] See [http://www.comsol.com] for details. 

[64] H. Busch, Annalen der Physik 64, 401 (1921). 

[65] Generally, to find a given solution and to trace it out 
over some current range, the current was first ramped 
from zero to high currents with a homogeneous temper- 
ature distribution. In order to get solutions with domain 
formation, a proper initial condition for the temperature 
profile of the ETD solution was made and, if successful, 
the current was increased or decreased in small steps, 
resulting in a branch with domain formation. 

[66] In other systems, stable configurations with few or even 
many hot spots in a sequence were considered for the 
case, when the out-of-plane resistivity is much higher 
than the in-plane resistivity [T] 1711 172] . In this case the 
high out-of-plane resistivity results in an effective self- 
consistent shunt that is stabilizing the hot spot structure. 

[67] K. L. Chopra, L. C. Bobb, and M. H. Francombe, Journal 
of Applied Physics 34, 1699 (1963). 

[68] See [http://cryogenics.nist.gov] for details. 

[69] C. Guerlich, S. Scharinger, M. Weides, H. Kohlstedt, 
R. G. Mints, E. Goldobin, D. Koelle, and R. Kleiner, 
Phys. Rev. B 81, 094502 (2010). 

[70] R. Werner, M. Weiler, A. Y. Petrov, B. A. David- 
son, R. Gross, R. Kleiner, S. T. B. Goennenwein, and 
D. Koelle, Appl. Phys. Lett. 99, 182513 (2011). 

[71] A. A. Akhmetov and R. G. Mints, Journal of Physics D: 
Applied Physics 16, 2505 (1983). 

[72] A. A. Akhmetov and R. G. Mints, Journal of Physics D: 
Applied Physics 18, 925 (1985). 



